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Abstract: We will review some known exact solutions for the steady state of the open quantum Heisenberg XXZ spin 
chain coupled to a pair of baths [1], The dynamics is modelled by the Lindblad master equation. We also review how to 
calculate some relevant physical observables and provide the statistics of spin current assuming the spin chain is weakly 
coupled to the baths [2]. 


1. INTRODUCTION 

In out-of-equilibrium physics the properties of the 
long time limit of systems are of paramount importance. 
Generically, systems when driven out-of-equilibrium re¬ 
turn to equilibrium after sufficiently long times. How¬ 
ever, cases when they do not are more interesting. There 
exists various setups which allow this. Here we will focus 
on systems driven out-of-equilibrium by a pair of infinite 
baths. The dynamics of the entire system+baths follows 
the standard Liouville-von Neumann equation. 


d_ 
d t 


p{t) = if p(t ) 


— i[-^sys+env, p{t )], 


(1) 


which is impossible to solve. However, since we are only 
interested in the dynamics of the system with Hamilto¬ 
nian, H, we can trace out the bath part. Then, if we in¬ 
voke certain approximations, namely the Born-Markov 
and rotating wave approximation and assume that the 
coupling between the system and the baths is weak, we 
are left with the Lindblad master equation for the reduced 
system dynamics [3], 


^Psys — if Psys — i 11 . p S ys] T Z ( Psys ) ■ 

^f(Psys) • ^ ' U; ^LkPsysLfo — ; Psys J'^ ■ 


( 2 ) 


where Z is the non-unitary part of the dynamics repre¬ 
senting the action of the baths on the system with so- 
called Lindblad jump operators Lj., and I i are the cou¬ 
pling constants. It will be useful to define (a,&H)p = 
[H, p], where ad H is a linear superoperator acting on the 
space of operators). Of particular importance is the fixed 
point, or steady state, p^. 


if Poo =-i ad iFpoo + #(poo) = 0. (3) 


We will focus on this state for the Heisenberg XXZ spin 
chain of size n given by the Hamiltonian, 

n 

n KXZ = X>5?i = 

J=i 

n 

= X 2 (°j + a J+1 + a 7 °i+i) + Aa^ + i,(4) 

j=i 

where er“ are the standard Pauli matrices acting on site j. 
Namely, 

tr“ = 1 2 j-l <8> cr“ (g) 1 2»-j, (5) 

and 1 k is the identity matrix of dimension k. The pa¬ 
rameter A is the anisotropy. 


2. QUANTUM INTEGRABILITY 


Referring to [4] consider the U q { sl 2 ) algebra, which 
is a quantum deformation of the 3(2 algebra and is de¬ 
fined by the standard q— deformed algebraic relations of 
its generators, s^, s+, sj, 

[ S X S 7 ] = [ 2s s] 9 > K, s t] = ±s t, (6) 


given by a Verma module representation, & s 

OO 

s* =X( s-fc )l fc ) ( fc l> 

k—0 

00 

k—0 

00 

=^2[2s- k] q \k + l) (k\, (7) 

k—0 


and characterised by a continuous spin parameter s, 
where [ x\ q := ( q x — q~ x )/{q — <Z -1 )- This then al¬ 
lows to define a general f7 g (sl 2 ) invariant Lax operator, 

L(<p, s ) G End (6 1/2 ® 6 S ), 


(sin (<P + 7 s s) (sin7)s s \ 
V (sin 7)3+ sin {tp — 7 s^) J 


( 8 ) 


where <p G C is the spectral parameter. Note that ©i/ 2 
corresponds to the physical space and 6 S to the “auxil¬ 
iary space”. Defining A = cos( 7 ), q = exp(i 7 ), and, like 



Eq. (4), ft xxz = 2 (<t + 0 p cr~ +cr~ 0 p a + ) + Act 2 0 p ct z 
( where 0 p refers to the tensor product in physical space), 
one can show that the Sutherland equation holds, 

[/i xxz , L(</>, s) 0 P L(<p, s)] = (9) 

2sin7(L(v?,s) 0 p L v (<p, s) - L p (</>,s) 0 p L(/?,s)), 

with 

s) := d v L(ip, s) (10) 

= cos ip cos ( 7 S Z ) 0 cr° — sin <p sin ( 7 S Z ) 0 er z , 

which is just the differential (in the spectral parameter ip) 
form of the famous Yang-Baxter equation. 

3. MAXIMALLY DRIVEN OPEN XX Z 
SPIN CHAIN 

Now we will present the results of Ref. [1] using a 
more compact form developed e.g. in Ref. [4], Recall the 
Lindblad master equation, Eq. (2) and take two Lindblad 
jump operator acting on boundary sites 1 and n of the 
spin chain, L \ and L n , representing a pure source and 
pure sink on the two ends, 

L 1 = y /eat, L n = y/ea~, (11) 

where e is the spin bath coupling. We are interested in the 
steady state of the system defined by Eq. (3), for which 
we take a now ubiquitous Cholesky-like ansatz, 

Poo = ss\ (12) 

with 


Now that we know the matrix product state form of the 
steady state explicitly, these observables can be com¬ 
puted efficiently numerically for arbitrary system size n, 
as well as, in certain cases, their exact analytical expres¬ 
sion in the thermodynamic limit, n —> oo. 

First let us define, following Eq. (12), a double Lax oper¬ 
ator L € End ( 0 i /2 0 & s 0 6 -s), 

L = Lg) a L t , (18) 

where 0 a refers to a tensor product of a pair of auxiliary 
spaces. 

For simplicity we will focus on only the two of the 
most relevant observables, the spin profile (cr|) and spin 
current, (j z ) = i((cr+a fc+1 - cr^a^ +1 )). To calculate 
these expectation values we will utilise three operators 
defined using Eq. (12), 

T = tr p L (19) 

V = tr p (<j“L) (20) 

W = tr p (i(cr + 0 p oj : +1 - o~ 0 p cr + )L 0 p L), (21) 

where the (partial) trace and the outer product are taken 
with respect to the physical space only. We will write 
|0> 0 a |0) := |0» Furthermore, we should define the trace 
of the un-normalised steady-state density operator as the 
non equilibrium partition function, 

= trpoo = «0|T"|0». (22) 


S = (0| L(<p, s)® n |0), (13) 

where L(<p, s ) is the Lax operator defined by Eq. ( 8 ). 
Plugging Eq. (12) into the equation for the steady state 
Eq. (3) and using the Sutherland equation Eq. (9) we are 
left with 10 boundary equations which we can solve for 
the spin parameter s and the spectral parameter <p. The 
solution is unique and is given by 

e = 4i[s ]~ 1 cos ( 7 s), (14) 

<p = 0. (15) 


The first useful relation between these operators, which 
can be shown using the algebraic relations Eq. ( 6 ), and 
which manifests the continuity equation for is, 

W = -2i[s] q T, (23) 

which leads to a result first seen in related classical asym¬ 
metric exclusion process models [5], 

(J) := (f k ) = - 2 i[s] q ■ ( 24 ) 


In the isotropic case A = 1 (7 —» 0) this reduces to, 

s = 4ie _1 . (16) 

Note that since s is purely imaginary (L(<p, s)p = 
(L (ip, -s)) T , (L (ip, -s)) T G End ( 6 1/2 0 6 _ s ) From 
now on, since we have calculated ip and s we will write 
for simplicity L := L (<p, s). Since we know the steady 
state density operator explicitly we can compute the ex¬ 
pectation values of observables in the long time limit. 

4. EXPECTATION VALUES OF 
OBSERVABLES 


We are interested in the expectation value of a local 
observable O in the long time limit. 


, Q s = tr(OPoo) 
trpoo 


(17) 


The operators T,V have effective finite rank m + 1 
for a dense set of anisotropies given by rational num¬ 
bers l/m, 7 = 7 rZ/rn, (recall that A = cos( 7 )) which 
densely cover the easy-plane regime |A| < 1. This al¬ 
lows for closed form expressions for small m where we 
can diagonalize these operators analytically. For exam¬ 
ple, for A = 1/2 = cos7t/3, in the thermodynamic 
limit, one finds flat spin profiles and can prove ballistic 
transport by explicitly computing the limit (see Fig. 1) 

/JA| ( v / 81+74e 2 +9e 4 -T-3e 2 )e 

\* 7 / |n—voo — 4(i_|_ e 2) • 

In contrast, in the the easy-axis regime | A| > 1, where 
the before-used effective truncation of the auxiliary space 
is not possible, one finds an asymptotically exponentially 
decaying current and hence insulating behaviour, (J) oc 
(| A| + -\/A 2 — \)~ n (Fig. lb-inset) with consistent kink¬ 
shaped magnetization profiles (Fig. la). 
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Fig. 1 (As taken from [1].) Spin profiles (a 7 j) at n = 100 

(a) , and spin currents (J) vs. size n (b), for A = 3/2 
(dashed), A = 1 (dotted/blue), A = 1/2 (full curves), 
all for three different couplings e = 1,1/5,1/25 us¬ 
ing thick, medium, thin curves, respectively. Red 
full curves show closed-form asymptotic results [see 
text]: (a 7 ) = cos tt , ( J) = n 2 e~ 1 n ~ 2 for A = 1 
in the main panels (a,b), and (J) ex e - narcoshA i n 

(b) -inset. 




Spin current is the flow of magnetization, M := 
^™_ 1 a 7 . That is, if we label the amount M transported 
in time t from the first bath to the second bath by N(t) 
then in the long time limit, J = lim^oo ^ N(t ). 

Magnetization is conserved by the Hamiltonian, but is 
changed by the Lindblad jump operators, Eq. (2). There¬ 
fore, the Lindblad jump operators are the ones which 
drive the current through the system. One approach to 
this problem, using master equations, which has attracted 
a lot of attention quite recently [6-11], is the method of 
so-called “full counting statistics” by utilising a “count¬ 
ing field”. In order to use this method first let us split the 
dissipator in the Lindblad master equation Eq. (2) into a 
jump part and a dissipative part, 

& UmP P ■■= E L mP L L 2> AiSS p ■= \ Ei4^ m , Ph 

fi m 

(28) 


Linally let us examine the critical isotropic point A = 
1, s = —. We will utilise a useful algebraic relation be¬ 
tween T and V, 


Depending on whether L m changes magnetization M by 
±1, i.e., [M,L m ] = ±L m , one modifies the jump part of 
the dissipator £A |u,np by introducing a counting field, x> 


[T, [T, V]] + 2{T, V} - 8s 2 V = 0. 
We will also need, 

2f n _ 2 / ( 4 n - 3) 2 


6Y 1 

°£n— ] 


~ £ 


32 tt 2 


— a + 1 + @(n ) 


(25) 


(26) 


where a is a constant, and the boundary conditions, 

«0|(T — V) = «0|, (T + V)|0)) = |0)). (27) 


®T P P - E e ±iX L m pLl (29) 

m 


which counts how many times a certain Lindblad jump 
operator has driven the current in either the left (+) or 
right (-) direction. With this modification in place one 
can show that the cumulants of the current are given by. 


<J m ) c :=fim ±([N(t)] m ) c 


d m \(x) 



(30) 


Then we can multiply the algebraic relation Eq. (25) by 
((0|TP _1 from the left, and T 71- - 7-2 ^}) from the right 
and use Eq. (26) and the boundary equations from Eq. 
(27) and in the continuum limit, M{x = 7ri_) := 
(cr'fj, find a second-order differential equation M"(x) - 
—7 t 2 M(x) + and solve with M(x) = cos 7 tx, or 

(crp ~ cos7r/5y- Also, Eq.(26) implies anomalous sub- 
diffusive scaling (J) « 7r 2 £ -1 n~ 2 to leading order in the 
thermodynamic limit. 

5. FULL COUNTING STATISTICS 

In this section we present the results from Ref. [2]. An¬ 
other interesting question concerning observables, with 
them being probabilistic in nature, is what are, not only 
their expectation values, but all other higher moments 
as well? The most relevant physical observable out-of- 
equilibrium is the current of some quantity. It is very 
challenging to study these fluctuation properties of the 
current which in general depend not only on the (asymp¬ 
totic) steady state of the system, but also on the correla¬ 
tions at earlier times. Therefore, our knowledge of the 
steady state density operator is not enough to compute 
the higher moments of spin current, J, even though it 
is enough to compute just (J), as shown in the previous 
section. 


Here A(x) is a leading eigenvalue (of maximal real part) 
of the modified Liouvillean, 


-iadiT + e(#f mp 



p(x) = Hx)p(x) 


(31) 


The full proof of this statement is rather long [12], how¬ 
ever it may be intuitively understood by observing a re¬ 
duced density matrix pAr(f), that is p(t) projected to a 
subspace of N spin-transfers between the two baths in 
time t. The trace of this, Pat(() = trpiv(t), is the 
probability of N spin transfers in time t. By perform¬ 
ing a Fourier transform (in N) of this reduced den¬ 
sity matrix, p(x, t) = PN(t)e~ lxN , it may then be 
shown (by observing the action of the generator of the 
time evolution) that the Lindblad master equation, Eq. 
(2), has the jump superoperator modified, ^ unlp p := 
E^e^L^pLl,,, so that it depends on the counting 
field x an d in which direction (//;/) a specific Lindblad 
operator L Mj „ drives the flow. Furthermore, if we nor¬ 
malize tr p(x, t = 0) = 1, the largest eigenvalue of the 
Liouvillian, A(x), corresponds to the cumulant generat¬ 
ing function for the current distribution in the long time 
limit [12], since for large t, p(x, t ) « e x ^ x ' )t p(x, t = 0). 
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Fig. 2 (As taken from [2]) An illustrative example of a 
parity symmetric spin system. The system is coupled 
to the two baths (represented by boxes), at site 1 and 
site n. The baths act via decoherent jump operators 
with corresponding rates (a, b, c, d). The arrows be¬ 
neath the jump rates indicate in which direction the 
spin current is being driven. 

When we take H = ll xxz \ Eq. (4) and maximal driv¬ 
ing driving, Eq. (11) finding A(x) analytically seems im¬ 
possible. However, we are in good shape since one of the 
assumptions when deriving the Lindblad master equation 
is that the system-bath coupling e is weak. This allows 
us to perform a perturbation series expansion in e. Quite 
remarkably, when we do this we find universal results, in 
the leading order for the full counting statistics of spin 
current, not just for H xxz , but for any parity symmetric 
Hamiltonian, H. This condition means that there exists 
an operator P, with P 2 = 1, such that PH = HP, and 
satisfying at least one of the following properties 

P<=<P, or Pal n = -al n P, (32) 

with an additional weak requirement, namely that also all 
conserved operators, Qk (ad II Qk = 0), are parity sym¬ 
metric [P, Qk] = 0. We can also drop the requirement of 
maximal driving via Lindblad jump operators, Eq. (11) 
and instead use the most general set of Lindblad jump 
operators acting on site 1 and n, 

I/ +;+ = y/eaa// , P+,- = Vebcr/, 

L_ t+ = y/sca+, L__=Vsda~. (33) 

See Fig. 3 for a graphic illustration. We will assume that 
the conditions of Evans theorem hold (or the Liouvillean 
can be symmetry reduced [13]) so the fixed point p{\) is 
unique. We perform a perturbation series in £ (using the 
fact that A(x) is an odd function of £), 

OO OO 

p(x) = ^2 (^) p p (p) (x), Hx) = '^£ 2p ~ 1 ^ {2p ~ 1) (x), 

p —0 p—1 

(34) 

Plugging this into Eq. (31) (defining S> x := ^ ump + 
3> diss ) one arrives to two equations for the first two orders 
in s, 

(adP> (0) = 0, (35) 

(ad fT)p (1) + J> x/ o (0) = A ( V 0) . (36) 


For Eq. (35) we take as an ansatz, 

n 

p(°) = 2-" XJ(1 +KxK). (37) 

j=i 

with free parameter i/(x), which may be found by requir¬ 
ing that — A^y 0 ) £ Im(adfT) (that is, there 

exists a solution to the first order equation Eq. (36)), as¬ 
suming the aforementioned parity symmetry requirement 
holds, Eq. (32). By taking the trace of the first order 
equation (36) we also find A*- 1 ), 

A (1) = 2a/ 1 + (e -2ix - 1 )ad + (e 2i * - l)bc - 2, (38) 
AW - (e~ ix - 1 )(a + d) + ( e ix - 1)(6 + c) 
(e~A — l)(a — d) + ( e ix — 1) (6 — c) 

We then apply Eq. (38) to calculate all the cumulants for 
this wide class of spin systems via Eq. (30). For instance, 
the expectation value of the spin current is = 

| (ad — be). Closed form expressions for higher cumu¬ 
lants were obtained in the same way, but are lengthy and 
therefore we will not write them out. However, they sig¬ 
nificantly simplify if we consider a symmetric driving in¬ 
stead of a general one (33), 

a = d={ I+aO/2, b = c = (1 - n)/2, (39) 

where the driving strength parameter // controls the non¬ 
equilibrium forcing due to unequal average spin polar¬ 
izations of the two baths. Then we have v = 0, = 

2“" 1, so A^ 1 ) = —1 + cos x — ift sin x, and (/ 2 ^ +1 ) c = 
Sfj,/2 for odd cumulants and (//() c = e/2 for even cu¬ 
mulants. Extreme driving p = 1 hence results in the 
Poisson distribution ( I m ) c = const. 

Under symmetric driving, Eq. (39), we may also find 
the third order correction to the current cumulant gener¬ 
ating function A 1 ''A for the XXZ spin chain II = // xxz . 
In order to do this we will turn again to the known so¬ 
lution for maximum driving without the counting held, 
Poo = SS\ where S is defined by Eq. (13). 

Define an operator Z from S by (see [14]) 

9 S 5| S= 0 = 27Sm7 ^ + 7COS (40) 

smy> 

where, as before A = cos( 7 ), ip is the spectral parameter, 
s is the continuous spin parameter and M is the magne¬ 
tization. Calculating d s ad iT xxz 5'| s= o using the Suther¬ 
land equation Eq. (9) one finds for Z, 


[H xxz ,Z)=al~a z n . (41) 

Remembering that for symmetric driving, Eq. (39) 
P( 0 ) = 2“” 1 we can immediately solve the first order 
equation Eq. (36), 2 n = c^(Z — Z^) (in a simi¬ 
lar way as when there is no counting held X = 0 [15], 
but multiplied by a different constant), c - 11 = (—p — 
p cos x + i sin x) /2. The second order equation, 

(ad if)/ 2 ) +i> x p (1) = A (1) p (1) - (42) 
















may be solved through a longer computation by taking an 
ansatz, similar in form to the one for \ = 0 [15], namely, 

2V 2) = c (1) 4 2) (Z - Zt) 2 - c (1) 4 2) [Z, Z f ], (43) 

and finding that c) ; = |(— /j — pcosx + i sin x) and 
cf = 1 (cos x — i/rsinx). The solution to the sec¬ 
ond order equation determines (r 2> only up to an ad¬ 
dition of a linear combination of conserved quantities 
Qk, ad 1 l xy: ' v (}k, he., the full solution can be written as 
p( 2 ) = p( 2 ) q- akQk- These turn out to be irrelevant 
however, because by taking the trace of the third order 
equation, 

A< 3 > = tr(4p^' - A(V 2) '). (44) 

and noting that the only relevant component of p <2> , 
£> x (crf — tr^), which contributes to the trace, is not a 
conserved quantity of ad // xxz due to the existence of 
a solution Z to Eq. (41). 

We therefore have, 

A (3) = (—+ cos x + isinx)tr[(—<7f + <T^)p (2) ]. (45) 


This finally allows us to write the third-order correct to 
the current cumulants. 


</(3)} c =-e 3 /(n) (9 ^ gg +1) > ( 4 6) 

2fc+l\ _ .3^„ ^(9 t+1 -l+3(9 fc -l)M 2 ) 




256(2fe+l)! 


where f(n) = (L| T" |R) — (L| T” -1 |R) and T is ex¬ 
actly the transfer matrix from Ref. [15], acting on auxil¬ 
iary Hilbert space with two ground-state vectors |R) and 
|L) (for details of how to derive this transfer matrix see 
also [14]). We can again compute f(n) efficiently for any 
anisotropy A of the form A = cos(7 rf/m), l,m £ Z. For 
instance, for A = 1, we have, f(n) = n— 1, and for A = 
1/2, f(n) = (5(—8)™ - 6(—5) n + 10). 

Finally we contrast our results with a spin system which 
does not fulfil the parity symmetry requirement from Eq. 
(32), namely a XXZ spin chain in a staggered field as 
shown in Fig. 3 


6. CONCLUSIONS 

We have shown how to find the steady state analyt¬ 
ically for the Heisenberg XXZ spin chain undergoing 
maximum driving under the Findblad master equation 
and compute the spin current, (J) and magnetization pro¬ 
files {u%). Depending on the value of the anisotropy A 
we find, in the thermodynamic limit n —> oo, either bal¬ 
listic (for a dense subset of |A| < 1), insulating (for 
|A| > 1), or anomalous transport (A = 1) with the spin 
current J decaying as J oc 1/n 2 . We have also shown 
what the universal results in leading order of perturbation 
in system-bath coupling e for the full counting statistics 
of spin current for any driving acting on the boundary 
sites 1 and n only are, as well as the next-to-leading order 
corrections for the Heisenberg XXZ spin chain undergo¬ 
ing symmetric driving. 



m 

Fig. 3 (As taken from [2]) The first four current cumu¬ 
lants obtained numerically for the XXZ spin chain 
and staggered field XXZ spin model, H = Hxxz + 
Sj=i h(±iy UjZ, with field strength h , for n = 4, 
A = 0.5, /r = 0.5, £ = 0.1. Dashed (chained) lines 
indicate analytical results up to second (fourth) order 
in £. 
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